% This code computes
% 1) each individual's choice-specific value of buying from each firm
% 2) for each firm, the probability that individuals keep the firm's contract in
%    each 2nd-period state
function [v, pr_stay] = consumer_v_ftn(p1, r2, N, naic_xi, lapse, N_y, ...
    rho, crra, beta_c, K, state_pr_s, S, sigma_c2, eulerc, B1_c, B2_c, T1, ...
    Delta_xi, alpha, D,...
    no_insurer_fe, ...
    n_variety, variety_u, delta, use_log_n, correct_beliefs,...
    resource_y1, resource_y2)

% predefine outcome variables
v = zeros(N_y, N.ist);
pr_stay = [];

% other data inputs
D_naic  = D.naic;
D_state = D.state;

% get data on additional agg risk
pr_mm  = D.pr_mm; % N.ist x M
eta_mm = D.eta_mm;
M = size(pr_mm, 2);

% loop thru each product jj
for jj=1:N.ist
    
    % jj's naic_code
    naic_jj = D_naic(jj);
    
    % jj's state (numeric value)
    state_jj = D_state(jj);
    
    % 2nd period probability distribution in state_jj
    pr_s = state_pr_s(state_jj,2:K+1)'; % K x 1
    
    % jj's prices
    p_jj1 = p1(jj, 1);
    r2_jj = r2(jj,:);
    
    % probability of additional agg states - depends on sales years
    pr_mm_jj = pr_mm(jj,:); % 1 x M
    
    % jj's unobserved characteristic
    % (1) persistent: at the insurer level: will be zero when no insurer FE.
    temp  = find(naic_xi(:,1)==naic_jj);
    xi_jj = naic_xi(temp, 2);
    
    % (2) transitory: at the insurer-state-year level
    Delta_xi_jj = Delta_xi(jj);
    
    % first-period annual consumption (N_y x 1)
    c1 = resource_y1-p_jj1;

    % jj's variety
    n_variety_jj = n_variety(jj);
    
    % second-period utility in each state k if lapse (N_y x K)
    U_lapse = S(state_jj).U0;
    
    % loop over mm==1,...,M
    if correct_beliefs==1
        
        EU_mat = zeros(N_y, M);
        
        for mm=1:M
            
            % 2nd-period prices conditional on mm
            p_jj2_mm = p_jj1*(1+0.01*r2_jj*eta_mm(mm));
            
            % second-period annual consumption in each state k if not lapse (N_y x K)
            c2 = repmat(resource_y2,1,K)-repmat(p_jj2_mm,N_y,1);

            % 11/2023: force termination if c2<1, i.e., cannot afford p2
            force_lapse = (c2<1);
            
            % second-period utility in each state k if not lapse (N_y x K)
            U_stay = B2_c*(alpha*u_c(c2,rho,crra) + Delta_xi_jj);
            
            % expected utility conditional on realized state (N_y x K)
            U = (1- repmat(delta, N_y, 1)).*U_stay + repmat(delta, N_y, 1).*U_lapse;
            U(force_lapse==1) = U_lapse(force_lapse==1);

            pr_stay(jj,mm).mat = 1-repmat(delta, N_y, 1);
            pr_stay(jj,mm).mat(force_lapse==1) = 0;

            % expected 2nd-period utility
            EU_mat(:,mm) = U*pr_s; %(N_y x 1)
            
        end

        % integrate over mm
        EU = sum(repmat(pr_mm_jj, N_y, 1).*EU_mat,2); % N_y x 1
        
    % if consumers believe stable rates, then mm won't affect utility
    elseif correct_beliefs==0
        
        % supply the initial rate as the 2nd-stage prices
        p_jj2 = repmat(p_jj1, 1, K);
        
        % second-period annual consumption in each state k if not lapse (N_y x K)
        c2 = repmat(resource_y2,1,K)-repmat(p_jj2,N_y,1);
        
        force_lapse = (c2<1);
        
        % second-period utility in each state k if not lapse (N_y x K)
        U_stay = B2_c*(alpha*u_c(c2,rho,crra) + Delta_xi_jj);

        % expected utility conditional on realized state (N_y x K)
        U = (1- repmat(delta, N_y, 1)).*U_stay + repmat(delta, N_y, 1).*U_lapse;
        U(force_lapse==1) = U_lapse(force_lapse==1);

        % expected 2nd-period utility
        EU = U*pr_s; %(N_y x 1)
        
        % not needed as computed in correct_beliefs=1         
        pr_stay = [];
        
    end

    v_jj = B1_c*(alpha*u_c(c1,rho,crra) + Delta_xi_jj) + (beta_c^T1)*EU;
    
    flag = c1;
    flag = (flag<1);
    flag = (sum(flag,2)>0);
    v_jj(flag==1)=-Inf; 

    % save
    v(:,jj) = v_jj;
end


